/**
 * 
 */
package ua.lviv.franko.solvers.old.twoDimention;

import java.util.ArrayList;
import java.util.Iterator;
import java.util.logging.Level;
import java.util.logging.Logger;

import org.jzy3d.maths.Range;

import ua.lviv.franko.integration.ConstantFunction2D;
import ua.lviv.franko.integration.Gauss2D;
import ua.lviv.franko.integration.IFunction2D;
import ua.lviv.franko.matrix.RowType;
import ua.lviv.franko.matrix.hash.lu.HashMatrixLUStyle;
import ua.lviv.franko.matrix.hash.lu.LUFactorization;
import ua.lviv.franko.mlsFunctions.oneDimension.weightFunctions.CubicSpline;
import ua.lviv.franko.mlsFunctions.twoDimension.BaseFunction;
import ua.lviv.franko.mlsFunctions.twoDimension.Particle;
import ua.lviv.franko.mlsFunctions.twoDimension.weightFunctions.IWeightFunction2D;
import ua.lviv.franko.mlsFunctions.twoDimension.weightFunctions.SquareSupport;
import ua.lviv.franko.solvers.twoDimention.EquationParameters;
import ua.lviv.franko.triangulation.SaveStructures.BaseFunctionValues;
import ua.lviv.franko.triangulation.SaveStructures.Node;
import ua.lviv.franko.triangulation.SaveStructures.Rectangle;
import ua.lviv.franko.triangulation.SaveStructures.ValuesOnTriangle;
import ua.lviv.franko.visualisation.EFGResult2D;

/**
 * @author marko
 * 
 */
public class Solver3 {
	// general data
	protected int						quadraturePointsNumber;
	protected Gauss2D					gau;
	protected EquationParameters		parameters;
	protected HashMatrixLUStyle				mainMatrix;
	protected double[]					b;

	// data for background integration;
	protected ArrayList<Rectangle>		recList;
	// protected ArrayList<Node> nodeList;
	// protected ArrayList<NodePair> edges;
	// data for base functions
	protected ArrayList<Particle>		particles;
	protected ArrayList<BaseFunction>	baseFunctions;

	public Solver3(int quadraturePointsNumber, ArrayList<Particle> particles, EquationParameters parameters, String folderDir, int PolinomialDegree, IWeightFunction2D weight) {
		// initialization
		this.quadraturePointsNumber = quadraturePointsNumber;
		this.gau = new Gauss2D(quadraturePointsNumber);
		this.particles = particles;
		this.parameters = parameters;

		// creation of base functions;
		this.createBaseFunctions(PolinomialDegree, weight);
		// continuation of initialization
		this.readData(folderDir);
		this.calculateValuesOnRects();
	}

	public EFGResult2D calculate() {
		for (int i = 0; i < b.length; i++) {
			for (int j = 0; j < b.length; j++) {
				System.out.print(mainMatrix.getValue(i, j) + "\t");
			}
			System.out.println();
		}
		for (int j = 0; j < b.length; j++) {
			System.out.println(b[j]);
		}
		for (int i = 0; i < b.length; i++) {
			BaseFunction bFunc = baseFunctions.get(i);
			if ((Math.abs(Math.abs(bFunc.getXi().getXi()) - 0.5) < 0.000001) || (Math.abs(Math.abs(bFunc.getXi().getYi()) - 0.5) < 0.000001)) {
				// if
				// ((Math.abs(Math.abs(bFunc.getXi().getXi())-0.5)<0.0000001)) {
				// for (int j = 0; j < mainMatrix.getN(); j++) {
				// mainMatrix.setValue(i, j, 0.0);
				// }
				mainMatrix.setValue(i, i, 100000.0);
				b[i] = 0.0;
			}
		}
		LUFactorization lu = new LUFactorization(mainMatrix);
		System.out.println();
		double[] res = lu.calculate(b);
		for (int i = 0; i < res.length; i++) {
			System.out.println(res[i]);
		}
		return new EFGResult2D(res, baseFunctions);
	}

	public void calculateValuesOnRects() {
		for (Iterator<Rectangle> iterator = recList.iterator(); iterator.hasNext();) {
			Rectangle rect = iterator.next();
			// get nodes in which to calculate;
			ArrayList<Node> tmpNodes = new ArrayList<Node>();
			ArrayList<BaseFunctionValues> bfvList = new ArrayList<BaseFunctionValues>();
			Node A = rect.A;
			Node C = rect.C;
			double ax = A.x, bx = C.x, ay = A.y, by = C.y;
			for (int j = 0; j < quadraturePointsNumber; j++) {
				for (int i = 0; i < quadraturePointsNumber; i++) {
					Node node = new Node((ax + bx) / 2.0 + (bx - ax) / 2.0 * gau.getX(i), (ay + by) / 2.0 + (by - ay) / 2.0 * gau.getX(j), false, 0);
					tmpNodes.add(node);
				}
			}
			// for (Iterator<BaseFunction> funcs = baseFunctions.iterator();
			// funcs.hasNext();) {
			for (int k = 0; k < baseFunctions.size(); k++) {
				BaseFunction function = baseFunctions.get(k);
				if (intersects(function.getXi(), new Particle((bx + ax) / 2.0, (by + ay) / 2.0, (by - ay) / 2.0))) {
					bfvList.add(null);
					continue;
				}
				int i = 0;
				BaseFunctionValues bfv = new BaseFunctionValues((quadraturePointsNumber) * (quadraturePointsNumber));
				for (Iterator<Node> iterator2 = tmpNodes.iterator(); iterator2.hasNext();) {
					Node node = iterator2.next();
					// double val = function.calculate(node.x, node.y);
					// double valX = function.calculateDerivatyX(node.x,
					// node.y);
					// double valY = function.calculateDerivatyY(node.x,
					// node.y);
					// bfv.setValue(i, val);
					// bfv.setDerivatyXValue(i, valX);
					// bfv.setDerivatyYValue(i, valY);
					bfv.setAll(i, function.calculateAll(node.x, node.y));
					i++;
				}
				// for(int j=0;j<quadraturePointsNumber;j++){
				// for (int i = 0; i < quadraturePointsNumber; i++) {
				// Node node = new Node((ax + bx) / 2.0 + (bx - ax) / 2.0 *
				// gau.getX(i),(ay + by) / 2.0 + (by - ay) / 2.0 *
				// gau.getX(j),false,0);
				// tmpNodes.add(node);
				// double val = function.calculate(node.x, node.y);
				// double valX = function.calculateDerivatyX(node.x, node.y);
				// double valY = function.calculateDerivatyY(node.x, node.y);
				// bfv.setValue(i, val);
				// bfv.setDerivatyXValue(i, valX);
				// bfv.setDerivatyYValue(i, valY);
				// }
				// }
				bfvList.add(bfv);
			}
			rect.values = new ValuesOnTriangle(tmpNodes, bfvList);
			// calculate function values and derivaties;
			//
			// for (Iterator<BaseFunction> funcs = baseFunctions.iterator();
			// funcs.hasNext();) {
			// BaseFunction function = funcs.next();
			// BaseFunctionValues bfv = new BaseFunctionValues(tmpNodes.size());
			// for (int i = 0; i<tmpNodes.size();i++) {
			// Node node = tmpNodes.get(i);
			// double val = function.calculate(node.x, node.y);
			// double valX = function.calculateDerivatyX(node.x, node.y);
			// double valY = function.calculateDerivatyY(node.x, node.y);
			// bfv.setValue(i, val);
			// bfv.setDerivatyXValue(i, valX);
			// bfv.setDerivatyYValue(i, valY);
			// }
			// bfvList.add(bfv);
			// }
			// triangle.setValues(new ValuesOnTriangle(tmpNodes,bfvList));
		}
	}

	public void generateMatrix() {
		mainMatrix = new HashMatrixLUStyle(RowType.ARRAY_TYPE, particles.size());
		b = new double[particles.size()];
		for (Iterator<Rectangle> iterator = recList.iterator(); iterator.hasNext();) {
			Rectangle rect = iterator.next();
			Node A = rect.A;
			Node C = rect.C;
			double ax = A.x, bx = C.x, ay = A.y, by = C.y;
			for (int i = 0; i < b.length; i++) {
				// System.err.println(i);
				if ((rect.values.getValues().get(i) == null)) {
					continue;
				}
				for (int j = 0; j < b.length; j++) {
					// if (intersects(i, j)) {
					// continue;
					// }
					// System.err.println(i+"\t"+j);
					if ((rect.values.getValues().get(j) == null)) {
						continue;
					}

					double tmp = gaussQuadratureIntegration(ax, bx, ay, by, rect.values.getValues().get(i).getFuncDerivatyXValue(), rect.values.getValues().get(j)
							.getFuncDerivatyXValue(), rect.values.getNodes(), parameters.fKxx)
							+ gaussQuadratureIntegration(ax, bx, ay, by, rect.values.getValues().get(i).getFuncDerivatyXValue(), rect.values.getValues().get(j)
									.getFuncDerivatyYValue(), rect.values.getNodes(), parameters.fKxy)
							+ gaussQuadratureIntegration(ax, bx, ay, by, rect.values.getValues().get(i).getFuncDerivatyYValue(), rect.values.getValues().get(j)
									.getFuncDerivatyYValue(), rect.values.getNodes(), parameters.fKyy)
							+ gaussQuadratureIntegration(ax, bx, ay, by, rect.values.getValues().get(i).getFuncDerivatyYValue(), rect.values.getValues().get(j)
									.getFuncDerivatyXValue(), rect.values.getNodes(), parameters.fKyx);
					tmp += gaussQuadratureIntegration(ax, bx, ay, by, rect.values.getValues().get(i).getFuncDerivatyXValue(), rect.values.getValues().get(j).getFuncValue(),
							rect.values.getNodes(), parameters.fPx)
							+ gaussQuadratureIntegration(ax, bx, ay, by, rect.values.getValues().get(i).getFuncDerivatyYValue(), rect.values.getValues().get(j).getFuncValue(),
									rect.values.getNodes(), parameters.fPy);
					tmp += gaussQuadratureIntegration(ax, bx, ay, by, rect.values.getValues().get(i).getFuncValue(), rect.values.getValues().get(j).getFuncValue(),
							rect.values.getNodes(), parameters.fQ);
					mainMatrix.setValue(i, j, mainMatrix.getValue(i, j) + tmp);
				}

				// double tmp =
				// gaussQuadratureIntegration(ax,bx,ay,by,rect.values.getValues().get(i).getFuncValue(),
				// rect.values.getNodes(), parameters.fF);
				b[i] += gaussQuadratureIntegration(ax, bx, ay, by, rect.values.getValues().get(i).getFuncValue(), rect.values.getNodes(), parameters.fF);
				// final int pos = i;
				//
				// for (int j = 0; j <
				// rect.values.getValues().get(i).getFuncValue().length; j++) {
				// System.out.println(rect.values.getValues().get(i).getFuncValue()[j]);
				// }
				// for (int j = 0; j < rect.values.getNodes().size(); j++) {
				// System.out.println(rect.values.getNodes().get(j));
				// }
				// System.out.println("-----------------");
				// double jaxz = gau.gau(ax, bx, ay, by, new IFunction2D() {
				//
				// @Override
				// public double calculateDerivatyY(double x, double y) {
				// return 0;
				// }
				//
				// @Override
				// public double calculateDerivatyX(double x, double y) {
				// return 0;
				// }
				//
				// @Override
				// public double calculate(double x, double y) {
				// double val = baseFunctions.get(pos).calculate(x,
				// y)*parameters.F(x, y);
				// System.out.println(baseFunctions.get(pos).calculate(x, y));
				// // System.out.println((x +"\t"+ y));
				// return val;
				// }
				// });
				// System.out.println();
			}
			System.out.println(rect.toString());
		}
		for (int i = 0; i < b.length; i++) {
			for (int j = 0; j < b.length; j++) {
				System.out.print(mainMatrix.getValue(i, j) + "\t");
			}
			System.out.println();
		}
		for (int j = 0; j < b.length; j++) {
			System.out.println(b[j]);
		}

	}

	private void createBaseFunctions(int PolinomialDegree, IWeightFunction2D weight) {
		this.baseFunctions = new ArrayList<BaseFunction>();
		for (Iterator<Particle> iterator = particles.iterator(); iterator.hasNext();) {
			Particle particle = iterator.next();
			BaseFunction bf = new BaseFunction(PolinomialDegree, particle, weight);
			bf.setNodeList(this.getParticleSupport(particle, weight));
			this.baseFunctions.add(bf);
			// bf.display();
		}
	}

	private double gaussQuadratureIntegration(double ax, double bx, double ay, double by, double[] phi1, ArrayList<Node> nodes, IFunction2D coefFunc) {
		double res = 0.0;
		int pos = 0;
		for (int j = 0; j < quadraturePointsNumber; j++) {
			double temp = 0.0;
			for (int i = 0; i < quadraturePointsNumber; i++) {
				Node node = nodes.get(pos);
				temp += gau.getW(j) * phi1[pos] * coefFunc.calculate(node.x, node.y);
				pos++;
			}
			res += temp * gau.getW(j);
		}
		res *= (bx - ax) * (by - ay) / 4.0;
		return res;
	}

	private double gaussQuadratureIntegration(double ax, double bx, double ay, double by, double[] phi1, double[] phi2, ArrayList<Node> nodes, IFunction2D coefFunc) {
		double res = 0.0;
		int pos = 0;
		for (int j = 0; j < quadraturePointsNumber; j++) {
			double temp = 0.0;
			for (int i = 0; i < quadraturePointsNumber; i++) {
				Node node = nodes.get(pos);
				temp += gau.getW(j) * phi1[pos] * phi2[pos] * coefFunc.calculate(node.x, node.y);
				pos++;
			}
			res += temp * gau.getW(j);
		}
		res *= (bx - ax) * (by - ay) / 4.0;
		return res;
	}

	private ArrayList<Particle> getParticleSupport(Particle p, IWeightFunction2D weight) {
		ArrayList<Particle> list = new ArrayList<Particle>();
		for (Iterator<Particle> iterator = particles.iterator(); iterator.hasNext();) {
			Particle particle = (Particle) iterator.next();
			if (weight.intersectedSupport(p, particle)) {
				list.add(particle);
			}
		}

		return list;
	}

	private void readData(String folderDir) {
		recList = new ArrayList<Rectangle>();
		int N = 5;
		double h = 1.0 / (double) N;
		double ax = -0.5;
		double ay = -0.5;
		for (int j = 0; j < N; j++) {
			for (int i = 0; i < N; i++) {
				Rectangle r = new Rectangle();
				r.A = new Node(ax + i * h, ay, false, 0);
				r.B = new Node(ax + (i + 1) * h, ay, false, 0);
				r.C = new Node(ax + (i + 1) * h, ay + h, false, 0);
				r.D = new Node(ax + i * h, ay + h, false, 0);
				recList.add(r);
			}
			ay += h;
		}

	}

	protected boolean intersects(int i, int j) {
		double Eps = 1E-7;
		Particle p1 = baseFunctions.get(i).getXi();
		Particle p2 = baseFunctions.get(j).getXi();
		boolean key = (Math.sqrt(p1.lengthTo(p2)) - (p1.getRho() + p2.getRho()) * Math.sqrt(2)) < Eps;

		if (!key) {
			Logger.getAnonymousLogger().log(Level.WARNING, "skipped");
		}
		return !key;
	}

	protected boolean intersects(Particle p1, Particle p2) {
		double Eps = 1E-7;
		boolean key = (Math.sqrt(p1.lengthTo(p2)) - (p1.getRho() + p2.getRho()) * Math.sqrt(2)) < Eps;

		if (!key) {
			Logger.getAnonymousLogger().log(Level.WARNING, "skipped in calculations of values");
		}
		return !key;
	}

	/**
	 * @param args
	 */
	public static void main(String[] args) {
		String FilePath = "C:/Users/marko/Documents/workspace/MyFEM/kvadrat/eps0";
		ArrayList<Particle> nodes = new ArrayList<Particle>();
		// параметри області:
		double ax = -0.5, bx = 0.5, ay = -0.5, by = 0.5;

		int n = 2, m = 2;

		double hx = (bx - ax) / (double) n, hy = (by - ay) / (double) m;

		// вузли:
		double rho = 1.5 * hx;
		for (int i = 0; i <= n; i++) {
			for (int j = 0; j <= m; j++) {
				Particle p = new Particle(ax + i * hx, ay + j * hy, rho);
				nodes.add(p);
			}
		}

		// equation parameters:
		EquationParameters params = new EquationParameters();
		params.setF(new ConstantFunction2D(10));
		params.setKxx(new ConstantFunction2D(1));
		params.setKyy(new ConstantFunction2D(1));
		params.setKxy(new ConstantFunction2D(0));
		params.setKyx(new ConstantFunction2D(0));
		params.setPx(new ConstantFunction2D(10));
		params.setPy(new ConstantFunction2D(10));
		params.setQ(new ConstantFunction2D(0));

		Solver3 solver = new Solver3(20, nodes, params, FilePath, 1, new SquareSupport(new CubicSpline()));
		solver.generateMatrix();
		solver.calculate().display(new Range(-0.5, 0.5), new Range(-0.5, 0.5));

	}

}
